This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

im <- readImage("papSmear.jpg")
dim(im)
[1] 493 335   3
plot(im) # raster method means within R

# reshape image into a data frame
df = data.frame(
  red = matrix(im[,,1], ncol=1),
  green = matrix(im[,,2], ncol=1),
  blue = matrix(im[,,3], ncol=1)
)
str(df)
'data.frame':   165155 obs. of  3 variables:
 $ red  : num  0.886 0.906 0.933 0.953 0.957 ...
 $ green: num  0.694 0.714 0.741 0.761 0.765 ...
 $ blue : num  0.71 0.729 0.757 0.776 0.78 ...
### compute the k-means clustering
K = kmeans(df,5)
str(K)
List of 9
 $ cluster     : int [1:165155] 1 1 3 3 3 3 3 3 3 3 ...
 $ centers     : num [1:5, 1:3] 0.885 0.502 0.966 0.707 0.61 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "1" "2" "3" "4" ...
  .. ..$ : chr [1:3] "red" "green" "blue"
 $ totss       : num 8016
 $ withinss    : num [1:5] 167 77.6 162.5 99.8 111.5
 $ tot.withinss: num 618
 $ betweenss   : num 7398
 $ size        : int [1:5] 18104 5536 105805 19219 16491
 $ iter        : int 6
 $ ifault      : int 0
 - attr(*, "class")= chr "kmeans"
head(K$centers)
        red     green      blue
1 0.8852673 0.6515596 0.7181804
2 0.5015577 0.3339716 0.4417304
3 0.9664161 0.8012112 0.8026724
4 0.7068195 0.6314813 0.6412725
5 0.6097929 0.5254669 0.5571331
df$label = K$cluster
head(df)
### Replace the color of each pixel in the image with the mean 
### R,G, and B values of the cluster in which the pixel resides:
# get the coloring
colors = data.frame(
  label = 1:nrow(K$centers), 
  R = K$centers[,"red"],
  G = K$centers[,"green"],
  B = K$centers[,"blue"]
)
dim(colors)
[1] 5 4
colors
# merge color codes on to df
# IMPORTANT: we must maintain the original order of the df after the merge!
df$order = 1:nrow(df)
df = merge(df, colors)
df = df[order(df$order),]
df$order = NULL
head(df)

Reshape our data frame back into an image

library(grid)
# get mean color channel values for each row of the df.
R = matrix(df$R, nrow=dim(im)[1])
G = matrix(df$G, nrow=dim(im)[1])
B = matrix(df$B, nrow=dim(im)[1])
  
# reconstitute the segmented image in the same shape as the input image
im.segmented = array(dim=dim(im))
im.segmented[,,1] = R
im.segmented[,,2] = G
im.segmented[,,3] = B
#im.segmented <- rotate(im.segmented, 90)
#im.segmented <- flop(im.segmented )
#dim(im.segmented)
#str(im.segmented)
#class(im.segmented)
#EBImage::writeImage(im.segmented, files = "im.segmented.jpg")
# View the result
#layout(t(1:2))
grid.raster(im.segmented)
#display(im.segmented, method = "raster", all = TRUE)
plot(im)

dim(im.segmented)
[1] 493 335   3
dim(im)
[1] 493 335   3
LS0tCnRpdGxlOiAiay1tZWFucyBjbHVzdGVyaW5nIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIGlzIGFuIFtSIE1hcmtkb3duXShodHRwOi8vcm1hcmtkb3duLnJzdHVkaW8uY29tKSBOb3RlYm9vay4gV2hlbiB5b3UgZXhlY3V0ZSBjb2RlIHdpdGhpbiB0aGUgbm90ZWJvb2ssIHRoZSByZXN1bHRzIGFwcGVhciBiZW5lYXRoIHRoZSBjb2RlLiAKClRyeSBleGVjdXRpbmcgdGhpcyBjaHVuayBieSBjbGlja2luZyB0aGUgKlJ1biogYnV0dG9uIHdpdGhpbiB0aGUgY2h1bmsgb3IgYnkgcGxhY2luZyB5b3VyIGN1cnNvciBpbnNpZGUgaXQgYW5kIHByZXNzaW5nICpDbWQrU2hpZnQrRW50ZXIqLiAKCmBgYHtyIGluY2x1ZGU9RkFMU0V9CiNzb3VyY2UoImh0dHA6Ly9iaW9jb25kdWN0b3Iub3JnL2Jpb2NMaXRlLlIiKQojYmlvY0xpdGUoIkVCSW1hZ2UiKQogbGlicmFyeSgiRUJJbWFnZSIpCmBgYAoKCmBgYHtyfQppbSA8LSByZWFkSW1hZ2UoInBhcFNtZWFyLmpwZyIpCmRpbShpbSkKcGxvdChpbSkgIyByYXN0ZXIgbWV0aG9kIG1lYW5zIHdpdGhpbiBSCmBgYAoKCgoKCmBgYHtyfQojIHJlc2hhcGUgaW1hZ2UgaW50byBhIGRhdGEgZnJhbWUKZGYgPSBkYXRhLmZyYW1lKAogIHJlZCA9IG1hdHJpeChpbVssLDFdLCBuY29sPTEpLAogIGdyZWVuID0gbWF0cml4KGltWywsMl0sIG5jb2w9MSksCiAgYmx1ZSA9IG1hdHJpeChpbVssLDNdLCBuY29sPTEpCikKc3RyKGRmKQpgYGAKCmBgYHtyfQojIyMgY29tcHV0ZSB0aGUgay1tZWFucyBjbHVzdGVyaW5nCksgPSBrbWVhbnMoZGYsNSkKc3RyKEspCmhlYWQoSyRjZW50ZXJzKQpgYGAKYGBge3J9CmRmJGxhYmVsID0gSyRjbHVzdGVyCmhlYWQoZGYpCmBgYAoKYGBge3J9CiMjIyBSZXBsYWNlIHRoZSBjb2xvciBvZiBlYWNoIHBpeGVsIGluIHRoZSBpbWFnZSB3aXRoIHRoZSBtZWFuIAojIyMgUixHLCBhbmQgQiB2YWx1ZXMgb2YgdGhlIGNsdXN0ZXIgaW4gd2hpY2ggdGhlIHBpeGVsIHJlc2lkZXM6CgojIGdldCB0aGUgY29sb3JpbmcKY29sb3JzID0gZGF0YS5mcmFtZSgKICBsYWJlbCA9IDE6bnJvdyhLJGNlbnRlcnMpLCAKICBSID0gSyRjZW50ZXJzWywicmVkIl0sCiAgRyA9IEskY2VudGVyc1ssImdyZWVuIl0sCiAgQiA9IEskY2VudGVyc1ssImJsdWUiXQopCmRpbShjb2xvcnMpCmNvbG9ycwpgYGAKCgpgYGB7cn0KCgoKIyBtZXJnZSBjb2xvciBjb2RlcyBvbiB0byBkZgojIElNUE9SVEFOVDogd2UgbXVzdCBtYWludGFpbiB0aGUgb3JpZ2luYWwgb3JkZXIgb2YgdGhlIGRmIGFmdGVyIHRoZSBtZXJnZSEKZGYkb3JkZXIgPSAxOm5yb3coZGYpCmRmID0gbWVyZ2UoZGYsIGNvbG9ycykKZGYgPSBkZltvcmRlcihkZiRvcmRlciksXQpkZiRvcmRlciA9IE5VTEwKCmhlYWQoZGYpCmBgYAoKIyMgIFJlc2hhcGUgb3VyIGRhdGEgZnJhbWUgYmFjayBpbnRvIGFuIGltYWdlCgpgYGB7cn0KbGlicmFyeShncmlkKQojIGdldCBtZWFuIGNvbG9yIGNoYW5uZWwgdmFsdWVzIGZvciBlYWNoIHJvdyBvZiB0aGUgZGYuClIgPSBtYXRyaXgoZGYkUiwgbnJvdz1kaW0oaW0pWzFdKQpHID0gbWF0cml4KGRmJEcsIG5yb3c9ZGltKGltKVsxXSkKQiA9IG1hdHJpeChkZiRCLCBucm93PWRpbShpbSlbMV0pCiAgCiMgcmVjb25zdGl0dXRlIHRoZSBzZWdtZW50ZWQgaW1hZ2UgaW4gdGhlIHNhbWUgc2hhcGUgYXMgdGhlIGlucHV0IGltYWdlCmltLnNlZ21lbnRlZCA9IGFycmF5KGRpbT1kaW0oaW0pKQppbS5zZWdtZW50ZWRbLCwxXSA9IFIKaW0uc2VnbWVudGVkWywsMl0gPSBHCmltLnNlZ21lbnRlZFssLDNdID0gQgojaW0uc2VnbWVudGVkIDwtIHJvdGF0ZShpbS5zZWdtZW50ZWQsIDkwKQojaW0uc2VnbWVudGVkIDwtIGZsb3AoaW0uc2VnbWVudGVkICkKI2RpbShpbS5zZWdtZW50ZWQpCiNzdHIoaW0uc2VnbWVudGVkKQojY2xhc3MoaW0uc2VnbWVudGVkKQojRUJJbWFnZTo6d3JpdGVJbWFnZShpbS5zZWdtZW50ZWQsIGZpbGVzID0gImltLnNlZ21lbnRlZC5qcGciKQojIFZpZXcgdGhlIHJlc3VsdAojbGF5b3V0KHQoMToyKSkKZ3JpZC5yYXN0ZXIoaW0uc2VnbWVudGVkKQojZGlzcGxheShpbS5zZWdtZW50ZWQsIG1ldGhvZCA9ICJyYXN0ZXIiLCBhbGwgPSBUUlVFKQpwbG90KGltKQpkaW0oaW0uc2VnbWVudGVkKQpkaW0oaW0pCmBgYAoK